Tuto2 • Tuto1 • Paper

Automated Discovery Tools to Study Behavioral Competencies of Biological Networks [📝 writing in progress]

Abstract

This paper presents a novel framework for mapping and exploring the behavioral competencies of gene regulatory networks. The framework leverages ideas from the behaviorist tradition to consider GRNs as entities acting in transcriptional space and provides a set of tools, leveraging computational models of curiosity-driven learning and exploration, for automated experimentation and behavioral characterization. We show that this framework is useful to discover the range of possible behavioral responses that the GRN can exhibit under different initial conditions, as well as to evaluate the robustness of the discovered behavioral abilities. After applying our framework on a database of curated GRN models, we discuss how the obtained "behavioral catalogs" can be a first stepping stone for better understanding the GRN functionalities as well as for designing drug-driven interventions in a biomedical or bioengineering context.

Introduction

As technology and biology continue to advance, we increasingly encounter a wide range of collective intelligences at every level of organization. However, unlike traditional agents such as humans and vertebrate animals, which have been extensively studied and characterized in terms of their behavior and responsiveness to external stimuli, there is a lack of a comprehensive "behavioral catalog" for these novel entities[1]. Recent work in basal cognition has shown that unconventional organisms, such as collectives of ants, bacteria, and slime molds, are capable of navigating complex environments to reach target states or locations despite encountering various obstacles in their way, akin to the navigation competencies observed in more traditional agents [21] [19] [25] [22]. But besides the navigation competencies observed in 3D space, Fields & Levin, 2022 [10] argue that it is also essential to start characterizing the competencies of agents operating in unconventional spaces such as for instance the transcriptional space (the space of possible gene expression levels), the physiological space (the space of possible physiological states and functions at the level of an organ or a body part), or even the morphological space (the space of possible morphologies of an organism). In fact, if we want to develop strategies for system-level biomedical interventions, we need to understand the goals and preferences of these unconventional organisms, how they operate and solve problems in their own space, and what stimuli we - as external observers - can use to drive the system towards desired outcomes, such as specific cellular functions or morphogenetic outcomes (Figure 1). For instance, by explicitly searching to discover diverse kinds of memory, Biswas et al. [4] [3] were able to reveal the behavioral abilities of biological GRN to store primitive forms of memories, including associative and transfer memories but also (undesired) forms of habituation and sensitization. Building on this knowledge, they were then able to find drug interventions to break those undesired behaviors [3].

In this paper, we propose a framework for studying the behavioral competencies of biomolecular networks, and in particular, of gene regulatory networks (GRNs). Gene regulatory networks are complex systems of genes and proteins that interact with each other, and that are responsible for turning on or off gene expression in response to internal and external stimuli. To ensure proper gene expression and cellular function, GRNs must solve a number of complex problems in the transcriptional space, including the regulation of gene expressions given changing environmental conditions or developmental cues, coping with noise and variability in gene expression, and ensuring robustness to perturbations and disturbances. To better understand the functionality of gene regulatory networks, significant effort has been made in their mathematical modeling. This has resulted in the development of large collections of publicly-available models, which were experimentally-determined and curated by biologists, before being shared on platforms such as the Biomodels website. Yet, even for these easily accessible and relatively simple numerical models, biologists often still lack of an effective understanding of the range of possible responses that the models can exhibit under different environmental conditions and different initial data.

In fact, building such a "behavioral catalog" is challenging for various reasons. First, the intervention space (the space of possible stimuli) is often very large and scientists do not know what interventions/stimuli in that space will be salient at the behavioral level, i.e. here will lead to changes at the transcriptional level. Moreover, because scientists do not know in advance the range of possible behaviors that GRNs can exhibit in their (latent) behavior space, coming up with strategies for (automatically) mapping and exploring the space of possible behaviors is not an easy task. Except for a subset of simple GRNs with only few and easy-to-find attractors, random screening techniques are inefficient due to the non-linearity and complexity of the dynamics involved. Machine learning (ML) techniques, on the other hand, are promising for assisting scientists in better mapping and navigating the high-dimensional parameter space of these biological networks [6] [31]. For example, recent work shown how evolutionary [12] [23] [28] and gradient-descent based [16] [27] optimization algorithms can help for the engineering of complex genetic circuits. It is noticeable however that most existing approaches follow two standpoints which, we believe, can restrict the scope of possibilities for applications in biomedicine and bioengineering. First, current approaches are mostly focused on rewiring the network structure (i.e. interventions at the hardware level) to achieve a desired functionality, typically for synthetic circuit engineering. In contrast, we believe that there is a lot of room for the development of stimuli based interventions at the software level, that could trigger (and exploit) the GRN problem-solving capacities (including its "latent" capacities that we might not be aware of), and that these could be applied both for synthetic engineered machines but also for in vivo biomedical applications. Secondly, these approaches generally consider the machine learning as an outcome-specific helper, typically using some optimization method toward a specific biosystem design goal. Not only those approaches are not straightforward to apply in practice (need a "proper" fitness function [12] and good-enough initialization to successfully navigate the chaotic optimization landscape) but, by focusing the search toward a specific (and often known to be possible) outcome they largely restrict the exploration of the solution space [27]. We believe that other (ML) tools can be more suited to first reveal and map the "latent" behavior space of these dynamical systems and then, possibly in combination with optimization-driven methods, exploit this knowledge to develop drug-driven interventions.

With this in mind, this paper proposes an empirical framework to conduct automated experimentation and behavioral characterization in ODE-models of biological networks (Figure 1).

First, for revealing the space of possible behaviors (pathways attractors) of the studied biological networks, we propose to use a family of curiosity-driven search methods, called intrinsically-motivated goal exploration processes (IMGEP) [2] [11]. Recent works have shown how IMGEPs can form useful discovery assistants in science, particularly for revealing behaviors of unfamiliar systems such as chemical oil-droplet systems [13], physical non-equilibrium systems [9] and models of continuous cellular automata systems [26][8][14]. Here, the curiosity-driven algorithm is used to automatically generate the sequence of interventions (GRN initial state) that will drive the exploration of the dynamical system by targeting a diversity of self-generated goals (Figure 1, column 3). Applying the method on a database of models from the BioModels website, we show that the space of possible behaviors revealed by the IMGEP algorithm (the discovered "behavioral catalog") is often much bigger that what we could have discover with simple intervention strategies like random screening.

Secondly, to go further than just mapping the space of possible behaviors of the dynamical system, we propose a variety of empirical tests to probe the "navigation competencies" of the studied biological networks. More particularly, we propose to assess the robustness of the discovered pathways to various perturbations such as noise and obstacle (Figure 1, column 4). This allows us to identify, among the discovered behaviors, various robust pathways consistently reaching "preferred" regions in transcriptional space despite encountering various perturbations on their way.

Altogether, we believe that they are several important avenues for reuses of such behavioral catalog (Figure 1, column 5). First, we discuss how it could provide useful insights on the degree of sophistication of the studied biological networks, an in particular on their abilities to have both robust and reprogrammable behaviors. Then, we discuss how it could help building a fuller picture of the energy landscape of the studied GRN [30] and how this could enable scientists to identify new (desired or undesired) pathways (as well as stimuli triggering those pathways). Lastly, we discuss how the acquired knowledge can serve as good initial basis for developing drug-driven interventions strategies.

As additional results, and to illustrate that our framework can be applied to other kinds of systems and problem spaces (as agnostic to their internal construction), we show how curiosity-search can be employed as an alternative (diversity-driven) strategy for the design of synthetic gene regulatory networks in a bioengineering context.

Finally, our codebase for simulation and experimentation on GRNs is written end-to-end in JAX, a high-performance python library which allows parallel experimentation and computational speedups, as well as automatic differentiation. We provide interactive tutorials and an executable version of the paper at https://developmentalsystems.org/curious-exploration-of-grn-competencies, as well as open-source python librarie SBMLtoODEjax (to convert SBML files into jax modules) and AutoDiscJax (to facilitate automated experimentation on these models).

Results

Analyzing GRN behavioral abilities in the transcriptional space: experimental setup

Table 1: Possible isomorphism between concepts in complex systems sciences and behavioral sciences.

Recently, there has been a growing interest in using tools and conceptualizations from dynamical systems in the behavioral sciences [15]. However, we believe that there could also be benefits to reverse this mapping and use concepts and tools from the behavioral sciences to characterize the behavior of dynamical systems. Indeed, we suggest that many dynamical systems, particularly those in developmental and synthetic biology, should be viewed as unconventional "organisms" that solve problems in their option space, rather than as passive collectives [1]. To illustrate this point, we propose an isomorphism between concepts in complexity and behavioral sciences in Table 1. This isomorphism raises scientific questions that are often overlooked when studying dynamical systems. For example, it explictly asks how we can identify the salient stimuli that these entities respond to, and what their goals and preferences are. More importantly, it allows for a reconsideration of our capacity for intervention, shifting our focus from micro-managing control to conditioning through stimuli, and making progress by constructing and sharing behavioral catalogs. Furthermore, this approach focuses on the functional analysis of behavioral capabilities and is independent of the organism's substrate or body, making it suitable for studying a wide range of subjects that solve problems in various spaces.

However, we also recognize that there are major roadblocks to applying traditional behaviorist approaches to unconventional areas such as developmental and synthetic biology. First, when dealing with unconventional organisms acting in unconventional spaces, scientists may not even know what the range of possible behaviors is. This limitation makes it challenging to design experiments to test the presence or absence of specific responses or preferences, and might limit the scope of discoveries. Second, behaviorism has traditionally relied on prior knowledge of salient stimuli, which allows them to test only a limited number of interventions or stimuli and compare results. Finally, behaviorism is traditionally depicted as never "opening the black box" but as the distinction between "internal" and "external" becomes blurrier, we might not want draw too strict distinction here.

In that prospect, we propose that the experimental setup should integrate more general definitions of the problem spaces (see Table 2 for a definition of those spaces and Figure 1, column 2 for a visualization of such spaces in the case of GRN). Experimentation should consider larger intervention spaces (including those that "open the black box") while prioritizing stimuli-driven interventions. We also suggest to recognize navigation in (unconventional) behavioral spaces as some form of problem-solving competencies and to acknowledge that we only know a small fraction of the reachable behavior space.

Table 2: Problem spaces

Here, as illustrated in Figure 2, our chosen experimental setup is the following. First, we chose the observation $o$ of the trajectory to store the gene expression levels of all nodes of the network model $(y_1(t), \dots, y_n(t))$ and for all simulation times $o=(\mathbf{y}(0), \dots, \mathbf{y}(T))$. Note that we can access and store those states as we use simulations, but we mainly do so for replay purposes as this choice does not really influences the discoveries (it is rather the choice of the behavior encoding $z=R(o)$ that does). Then, we chose the behavior descriptor $z$ to encode the last point of the system trajectory for two (target) nodes $z=(y_i(T), y_j(T))$. The encoded behavior represents the attractor (stable endpoint) that the trajectory is converging to in transcriptional space as, for long-enough simulation time $T$, the observed trajectories have typically converged to a fixed point (or reached a periodic orbit but we discard those cases, as explained in Materials and Methods). Here, for simplicity purposes, we only observe the final behavior of two nodes but obviously many other choices could be made for the behavioral embedding (e.g. ones that also characterize the dynamics of the rollout). Note that the extent of the behavior space $Z \in \mathbb{R}^2$ is unknown a priori, and depends on how much the GRN constrains the state space. Finally we chose the intervention to set the system's initial state $i=y(0)=(y_1(0), \dots, y_n(0))$ of the rollout trajectory. Once again, many other choices could be made (e.g. ones shown in Figure 2 d-e), but we wanted the intervention to be minimal in order to study the "native" responses of the studied networks (and influencing the network dynamical state instead of the network's parameters and/or structure). The extent of the intervention space $I \in \mathbb{R}^n$ is set by the experimenter, here we set it to $y_i(0) \in [r*min(o_{default}[i]), \frac{1}{r}*max(o_{default}[i])]$, where $o_{default}$ is the default trajectory (i.e. the one starting from the initial conditions given by the SBML file) and $r=20$ determines the scaling of the intervention space (see appendix X for an alysis of the impact of $r$, and for a proposed variant where the extent is not fixed a prior). Note that, because we have no prior knowledge of salient regions in the intervention space, we opted for such definition of the intervention space but this leads to relatively vast (and high-dimensional) spaces to explore.

Given such setup and an experimental budget $N$, we can make advantage of numerical simulations to explore the system responses under different interventions $i_1, \dots, i_N \sim I$ and construct a "behavioral catalog" $\mathcal{H} = \lbrace (i_k,o_k,z_k), k=1 \dots N \rbrace$ with the observed responses (see Materials and Methods for more details on how we perform numerical simulations) . One central challenge though, when dealing with high-dimensional intervention spaces, is that we do not know how to organize the sequence of interventions/stimulis nor how this will impact the discovered behavioral catalog. For example, as we can see in Figure 2 (c-f), large portions of the intervention space can lead to similar effects in behavioral space.

Automating the behavioral catalog construction with curiosity-driven search

To generate the sequence of experiments that will organize the exploration of the GRN (and hence the construction of the behavioral catalog), we propose to use a family of algorithmic processes called intrinsically-motivated (or curiosity-driven) goal exploration processes (IMGEP). Inspired by research in developmental psychology and neuroscience, and originally applied in the field of developmental robotics, IMGEP have shown how artificial agents (robots or AI) endowed with intrinsic motivation incentive (instead of external reward) are able to autonomously explore their environment and learn what effects can be produced by their actions [2][11]. More rencent applications have explored the use of IMGEPs as scientific "discovery assistant", where instead of controlling the (body) movements of a robotic agent the IMGEP was used to map the space of possible outcomes of complex self-organizing systems by generating the sequence of input parameters sent to the system [13][26][8][9][14]. To generate the sequence of actions/interventions that will drive the exploration of the environment/dynamical system, the IMGEP has several internal modules and key mechanisms that it uses to represent, self-generate and achieve goals based on based on intrinsic motivation signals (e.g. to promote diversity or learning progress) and previous experience (knowledge acquired throughout exploration).

Here, we use the IMGEP to iteratively select the GRN initial states (interventions $i$) that will condition the discoveries in $\mathcal{H}$, with the aim to maximize the diversity of discovered behaviors $z \in Z$ at the end of exploration, and given a limited budget of N experiments. Figure 1 (c) provides a pseudocode of the IMGEP pipeline, and we refer to Material and Methods for more details on the chosen IMGEP modules as well as to Part III of tutorial 1 for a step-by-step walkthrough of the method.

We then compare the discoveries made by our curiosity-driven algorithm (IMGEP) with the discoveries made by random search (RS), on a database of 432 systems collected on the BioModels platform [20]. For the IMGEP, we compare two variants one which has only m=1 optimization run per target goal and the other one which has m=3 trials per target goal (but the same overall budget $N$). Random search generates interventions $i_1, \dots, i_N \sim \mathcal{U}(I)$ by uniform sampling of the intervention space, given an experimental budget $N=900$ which is twice as much as the budget we give to the IMGEP variants (N=450). Our database of 432 system, where a system is a (model, observed nodes pair) tuple, have been created to discard systems that are "too easy" to explore such as ones that have very few easy-to-find attractors (see Materials and Methods for more details).

Figure 3 shows the resulting diversity discovered by the different algorithm variants, averaged over the 432 systems of the database. We measure diversity as the area of the discovered reachable space (see Material and Methods for more details on the diversity metric). We can see that the IMGEP variants are significantly better at revealing the diversity of possible behaviors in the explored systems when compared to random search, even when given twice as less budget. There is no significant difference however between the IMGEP variants m=1 and m=3. Figure 4 details the discoveries made by the different exploration algorithms for the random exploration strategy "worst cases" (systems (69,(0,4)), (455,(4,5)) and (647,(0,6))) and "best cases" (systems (272,(2,3)), (240,(1,4)) and (641,(1,2)). In the first 3 examples, we can see RS is really not efficient in discovering the reachable space in $Z$ whereas IMGEP is much better despite being given twice as less budget. The other 3 show some examples of "easy to explore" systems where RS is already quite efficient, but so is the IMGEP. Those results suggest that, when we do not know in advance the space of possible behaviors of the system, IMGEP is an interesting tool to efficiently (and automatically) assist the discovery of interventions in $I$ that lead to diverse effects in $Z$.

Moreover, IMGEPs are designed to be flexible and adaptable (see Appendix X), so they can be used in a wide range of applications and systems (see section X). We believe that such diversity-driven approaches are promising approaches to construct a first "catalog" of possible behaviors in unconventional entities, like GRN.

Assessing the robustness of the discovered behavioral abilities with empirical tests

For characterizing the "goals" of the GRN and the competency with which it can reach them, we designed a battery of tests to see wether the GRN were able to consistently reach the same endpoints despite encountering various perturbations in their way. More precisely, we tested how robust the GRN trajectories were to (i) noise to the gene level expressions, (ii) sudden perturbations or "push" in transcriptional space during its traversal toward endpoint, and (iii) energy barriers or "walls" that act as force fields in transcriptional space.

As detailed in Materials and Methods, those perturbations are implemented as parametrized distributions that we can sample from $u \sim U$ and that, together with the intervention $i$, influence the numerical simulation trajectory $o$. We have three family of perturbation, the noise perturbation $U_1$ that is parametrized by its standard-deviation (how much noise) and period (how often), the push perturbation $U_2$ that is parametrized by its magnitude (how big) and number (how many) and the

Our procedure, illustrated in Figure 1 - column d, is the following. We define 3 family of perturbations: the noise perturbation $U_n(\sigma_n, p_n)$ which is parametrized by its standard-deviation and period, the push perturbation $U_p(m_p, n_p)$ parametrized by its magnitude and number of occurrences, and the wall perturbation $U_w(l_w, n_w)$ parametrized by its length and number. For each family, we chosen 6 set of parameters, resulting in 18 distributions that we can sample from: $U_n(\sigma_n, p_n)$ with $(\sigma_n, p_n) \in \lbrace(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)\rbrace$, $U_p(m_p, n_p)$ with $(m_p, n_p) \in \lbrace(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)\rbrace$ and $U_w(l_w, n_w)$ with $(l_w, n_w) \in \lbrace(0.05,1),(0.1,1),(0.15,1),(1,0.1),(2,0.1),(3,0.1)\rbrace$.

Then, for each system $(I,Z)$ and its corresponding behavioral catalog $\mathcal{H}$, where $\mathcal{H}$ here is the one constructed by the imgep m=1 algorithm variant, we retrieve $K$ trajectories $\lbrace (i_k,o_k,z_k) \in \mathcal{H} \rbrace$. The $K$ trajectories are representatives discoveries, i.e. ones that cover well the reachable space, and we chose $K=45$ (10\% of N, as K=N would require a lot of compute). For each trajectory $(i,o,z)$ and for each perturbation distribution $U$ we sample $p$ perturbations $u \sim U$. For each perturbation $u$, we re-run the trajectory starting from the same initial state ($i$) but with the sampled perturbation $u$ acting on the trajectory. Then, we measure the distance between trajectory's endpoint prior and after perturbation, and average the results over the $p$ perturbations which gives us a proxy of what we call the sensitivity of the trajectory. Here, for each family we sampled $p=3$ perturbations which results in a total budget of $K*18*3$ rollouts per system. We refer to Materials and Methods for more details on the selection of the representative discoveries, on the robustness tests, and on the sensitivity metric.

Results are illustrated in Figure 4 where, for each of the 432 systems, we show the median sensitivity of the 45 representative discoveries over the different perturbation families. Overall, even though we do observe some degree of sensitivity to perturbations, the sensitivity remains relatively low at least with respect to what we could expect knowing that the systems we study here are not the traditionally-studied GRNs with only few strong attractors but systems that are in fact capable of many diverse behaviors, as shown in previous section. The figure also shows some of the "extreme" responses in term of robustness, i.e. some trajectories of the least and most robust systems to the tested perturbations. Very interestingly, this procedure allows to identify the systems that demonstrate quite remarkable "navigation competencies" and that, not only are able to robustly reach their endpoint despite several perturbations but show quite complex trajectories which is reminiscent of goal-reaching and delayed-gratification behavior that we observe in more conventional organisms.

Reusing the discovered behavioral catalog

These kinds of analysis can provide insights on 1) the latent functionalities and resilience abilities of the studied biological networks 2) on their degree of sophistication in comparison to one another and 3) our ability to intervene on the molecular pathways.

To have better picture on the GRN functionality

To develop insights on the degree of sophistication of the different GRNs

To develop therapeutic interventions

Adapting the framework to other systems and problem spaces

Discussion

Limitations:

• Models might be limited (and behaviors we discover might not translate in in vitro studies)
• Simulation uses jax ODE solver method but choice of the solver has influence on results [29]
• Here we use predefined behavioral space, which is limiting for more complex models.

Future work:

• Learning representations for behavioral characterization as in [26][8]
• Deploying similar framework but at the multi-cellular level [18]
• Applications in the lab like [24]

Materials and Methods

Numerical simulation

ODE (ordinary differential equation) models of gene regulatory networks represent each gene by a continuous variable (representing the gene's expression level), and model the interactions between genes as a system of differential equations, where the rate of change of each gene's expression level is a function of the expression levels of the other genes in the system. These ODE models are usually stored and exchanged using the Systems Biology Markup Language (SBML), a standard format for representing mathematical models of biological systems, including gene regulatory networks. SBML files contain information about the variables, parameters, and equations that describe the behavior of the system (as well as additional metadata about the model), are written in XML format, and can be found on online platforms such as the BioModels database. There is no single standard way to do numerical simulation of those models, as the XML files can be read and simulated by many software tools, and the choice of the simulator should be made based on one's needs and experience with programming.

In this work, we use a recently-developped python package called SBMLtoODEjax, which provides a set of tools to automatically parse and convert the SBML files into python modules. Those modules are then used to obtain time-dependent solutions for the expression levels of each gene, given a set of initial conditions and model parameters, by numerical solving of the differential equations. To facilitate automated experimentation and simulation of those python modules, we developed the AutoDiscJax python package, which provides a handful of other python modules and pipelines to intervene on the simulation dynamics (e.g. simulate genome and/or drug interventions), but also to automatically organize experimentation via different possible exploration approaches (such as random, optimization-driven and curiosity-driven search) and to perform various sorts of perturbations on the discovered trajectories (such as noise, pushes and walls on the grn states but also on the grn kinematic parameters).

We refer to tutorial 1 for a step-by-step walkthrough on the uses of these libraries to apply the proposed framework on an example GRN model [^hyunMathematicalModelingInfluence2003], as well as to tutorial 2 for another example use-case of the AutoDiscJax python package and proposed exploration algorithms for the discovery of synthetic gene circuits (presented in the Results section 5).

Database creation

We use the publicly available database of manually curated ODE models from the BioModels website [20]. Through the 1048 curated models that are referenced on the website, only 188 are usable for us for simulation due to diverse simulation errors or to too long compute times. As we are observing pair of target genes/nodes, our database actually consists of $\lbrace(model, (node1, node2))_i\rbrace_{i=1 \dots N}$ tuples, that we call system. We apply additional filters to filter out pair of nodes that xxx, which results in our final database of 432 systems. Appendix X provides additional details on the selection procedure and on the selected systems (biomodel identifiers, number of species, etc.).

Curiosity-driven exploration

The IMGEP internal modules include a goal-embedding module $R$ , which encodes the observation $o$ into the IMGEP goal space $Z$, a goal generator module which iteratively samples goals from the goal space $g \sim \mathcal{G}$, a goal-conditionned loss which defines a loss function that measures the agent's progress towards achieving its current goal, and a goal-conditioned optimizer which select promising intervention parameters $i$ and optimizes them (given a certain budget) based on the current goal and the goal-conditioned loss.

Note that we've made relatively simple choices for the implementation of the IMGEP goal representation (predefined encoding of $o$ into $z$), goal generation (uniform sampling of $g$ in the bouding hypercube of current set of points in $Z$) and goal-conditioned optimization (local diffusion and selection of $i$) strategies, but many other choices could be envisaged.

Robustness tests

Additional results